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Abstract. In an attempt to extend the range of model jamming transitions, we simulate systems of athermal 
particles which attract when slightly overlapping. Following from recent work on purely repulsive systems, 
dynamics are neglected and relaxation performed via a potential energy minimisation algorithm. Our 
central finding is of a transition to a low-density tensile solid which is sharp in the limit of infinite system 
size. The critical density depends on the range of the attractive regime in the pair-potential. Furthermore, 
solidity is shown to be related to the coordination number of the packing according to the approximate 
constraint-counting scheme known as Maxwell counting, although more corrections need to be considered 
than with the repulsive-only case, as explained. We finish by discussing how the numerical difficulties 
encountered in this work could be overcome in future studies. 

PACS. 45.70.-n Granular systems - 64.70.Pf Glass transitions - 82.70.Gg Gels and sols 



1 Introduction 

The mechanical properties of athermal aggregates such 
as granular media (e.g. sand), wet foams and emulsions 
are ultimately related to the elastic interactions between 
their respective constituents (grains, bubbles, droplets) 
and the morphology of the contact network [TJUKS] . How- 
ever, in contrast to chemical bonding in equilibrium sys- 
tems, where the network morphology can be derived [4], 
persistent elastic contacts form as a result of kinetic energy 
dissipation leading to arrest, resulting in a nonequilibrium 
state that is often referred to as 'jammed'. The mechan- 
ics of the arrested state is thus intrinsically linked to its 
history and loading, potentially allowing for unusual sus- 
ceptibility properties when different loads are applied [5]- 

One approach to understanding the nature of this pro- 
cess is to study model systems. Simulations by O'Hern et 
al. on spheres with Hertzian or truncated Hookean inter- 
actions jHlEI confirmed earlier predictions that the tran- 
sitional state delineating solid ('jammed') from non-solid 
systems is isostatic [8], this term referring to a tenuous 
solid where severing just one contact would initiate global 
break up. Macroscopic quantities such as pressure and 
elastic moduli were shown to scale above this transition 
with simple exponents that are independent of dimension. 
The interpretation of the transition remains unclear: a 
field-theoretic approach suggests a mixed nature [5], al- 
though a possible central role for elastic stability has also 
been suggested [lOlfTT] . 

Greater insight into jamming as a whole could be gained 
by expanding the range of transitions studied. This is the 
purpose of this paper: to investigate a different class of 



jamming transition using similar techniques to previous 
work. The key difference here is that the particle inter- 
actions are attractive when just touching, allowing for a 
transition to disordered solid dominated by tensile contact 
forces. Indeed, this is precisly what is found: a well-defined 
transition to a tensile solid at a lower volume fraction than 
the repulsive-only case. This is not an entirely abstract 
exercise, as non-Brownian particles will often have short- 
ranged attractions, such as from van der Waals [12] or 
surface tension [13] effects. Therefore probing the role of 
jamming in attractive systems will also help to understand 
the mechanics of a large class of materials. 

A visual demonstration of the effects of introducing at- 
tractive interactions is shown in Fig. [T] These snapshots, 
generated by the freely available demonstration package 
SandBox [14j . show polydisperse, 2-dimensional systems 
of circular particles with the same mixed attractive/re- 
pulsive potential and initial conditions as used in the rest 
of this paper, simulated via a dynamical rule similar to 
that of Kasahara and Nakanishi [15]. The message from 
this figure is that attractive systems tend to cluster into 
a heterogeneous network with large regions devoid of par- 
ticles, unlike the much more homogeneous repulsive-only 
system (also shown). It should therefore not be surprising 
that solidification occurs at lower densities. Repeated sim- 
ulations clearly demonstrate the system relaxes to a static 
state with a negative pressure for a range of densities. 

For simplicity, and to facilitate a greater correspon- 
dence with previous works, the results presented in the 
remainder of this paper were generated by an energy min- 
imization algorithm which has no proper mapping to a 
damped system (indeed, it has no time scale). It is there- 
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Fig. 1. (Top) Example of a relaxed 2-dimensional configura- 
tion of athermal particles at mechanical equilibrium, with a ra- 
dial interaction that is attractive for small overlaps (white lines, 
thickness proportional to log of force magnitude) and repulsive 
for large overlaps (black lines). For comparison a repulsive-only 
system at higher density is shown in the lower panel. 

fore important to stress that, although we find net-tensile 
states in the N — > oo limit with this algorithm, we as 
yet have no strong evidence this will remain true with 
dynamical relaxation rules. (The demonstration package 
described above is limited to systems of around N = 200 
particles) . Determining the thermodynamic limit with dy- 
namical relaxation rules is, however, beyond the scope of 
this article, and must await future verification. 

This paper is arranged as follows. In Sec. [2] the model 
and simulation method are introduced. The results are 
detailed in Sec. [3] in three parts: firstly the variation of 
pressure with volume fraction, which shows the existence 
of the tensile solid. Secondly the variation with system size 
around this transition point, to confirm the transition is 
sharply-defined in the infinite system limit. Finally the 
variation of various coordination numbers with density 
and how this relates to constraint counting measures of 
rigidity. Further discussion, including possible directions 
for future work, is provided in Sec. [4] 



2 Model 

An obvious starting point when modelling attractive par- 
ticles is to select a suitable interaction potential for the 
system in question. However, this runs the danger of being 
non-generic, and in the case of surface tension, the stan- 
dard potentials in use, such as JKR (Johnson, Kendall and 
Roberts) and DMT (Derjaguin, Muller and Toporov) [15] . 
are non-conservative and hence inappropriate for the en- 



ergy minimisation algorithm to be used. Therefore a sim- 
ple shifted and truncated Hookean potential will be used. 

The d = 3 dimensional system comprises of point par- 
ticles interacting via with a purely radial pair potential 
U(r), with r the scalar particle separation. U(r) is strictly 
zero for particles separated by distances r > 2R, with the 
maximum interaction range R is the same for each parti- 
cle. In terms of the overlap e = 2R — r, 

jjf \ _ f \k(e — £ s hift) 2 — U co : e > , s 
[) \ : e<0 [ > 

where the energy offset U co = ^ke^ hift to ensure continuity 
at e = 0. Thus, there is an attractive regime for small 
overlaps < e < £ s hift and a repulsive regime for larger 
overlaps e > £ s hift- A repulsive-only truncated Hookean 
interaction is recovered by setting £ s hift = 0. 

The system is prepared according to an easily repro- 
ducible T — oo quench protocol: N particles are added at 
random to an L x L x L cubic cell, where L is chosen to 
give the required volume fraction <j> = Nv p /L 3 with v p the 
volume of a single particle (note that this 'bare' volume 
fraction double-counts overlaps). For the data presented 
below, is calculated on the basis that the true radius of 
the particle is the minimum of U(r), but no changes in any 
trends were observed when using the maximum extent R. 

The total potential energy V({Xi}), with Xj the coordi- 
nates of the particle centres, is then minimised via a high- 
dimensional optimisation algorithm, either conjugate gra- 
dient or Broyden-Fletcher-Goldfarb-Shanno (BFGS) [17] . 
which were found to agree except near the rigidity tran- 
sition, for reasons to be discussed in Sect. [4] Although it 
may be tempting to associate energy minimisation with 
over-damped dynamics, it should be noted that these al- 
gorithms do not move along lines of steepest descent (which 
is computationally inefficient) . Thus strictly speaking there 
is no precise physical interpretation of this energy minimi- 
sation procedure. 



3 Results 

3.1 Variation of pressure with volume fraction 

For purely repulsive potentials, the pressure P in the final, 
relaxed state (which is straightforward to extract from the 
contact forces [IB]) is identically zero below some critical 
volume fraction (j) c ks 0.639, which is sharply defined in 
the N — » oo limit [BJ. No particles are touching in this 
state. For <fi > <f> c , however, the system relaxes to form a 
disordered solid with P > and finite elastic moduli. 

The introduction of attraction U co > changes this 
picture in two qualitative ways, as seen in Fig. [2j 

(i) P is zero below some critical density 4> c (U co ) that 
depends on U co , and is negative immediately above this. 

(ii) At some higher density (which also depends on 
Uco), P becomes positive. 

Numerical convergence was noticeably slower around <f> c , 
which is not surprising as the repulsive-only system is 
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Fig. 2. Plot of normalised pressure P/(k/R) versus the volume 
fraction 4> for U co = (repulsive-only), 10~ 3 (so £ s hift/27? « 
5%) and 1CT 2 (e s hift/27? « 15%). The system size was N = 
1000 particles. (Inset) A magnified region of the data, showing 
that the U co = 10 -3 data becomes negative. 



highly correlated here [19], and the attractive case is ex- 
pected to be also. Such correlations increase the algorith- 
mic demands in finding a local minimum. No such diffi- 
culties arose at the higher density where P crosses over 
from negative to positive, suggesting there are no unusual 
susceptibility properties associated with this P = point. 

It is interesting to note that the data for the narrower 
attractive well, U co — 10 -3 , attains net positive pressure 
P > at a lower volume fraction than the strictly repul- 
sive case U co — 0. This remains true (to a lesser extent) 
even when the maximum extent R of the particle is used 
the calculate the particle volumes, rather than the mini- 
mum of U (r). Now, given that the P < regime is weakly 
gelled (in the sense that \P\ is small), it is conceivable 
that when analysing noisier data such as that gained from 
experiments, the P < regime would be entirely missed 
and 'jamming' ascribed to the point when P becomes pos- 
itive, which is below random close packing 0r,cp ~ 0.64. 
Thus even a narrow attractive regime could significantly 
contribute to the phenomenon of random loose packing. 




Fig. 3. Fraction of 'jammed' systems (i.e. with versus 
4> for Uco = 10~ 2 and the system sizes shown. The solid lines 
are fits to 4{1 + tanh[(c/> — <^ m id)/a"]}. 




Fig. 4. m id — (poo and a for the fits in Fig. [3] plotted against 
1/N. The upper solid line is a fit to J3j|. No convincing fit was 
found for the noisier a data, including power law, exponential 
and logarithmic trial forms; the supplied line is a guide to the 
eye (in fact a fit to a cx ~N~ V with the same v as in (0). 



tern size TV, thus allowing the N — > oo limit to be extrap- 
olated. This is shown in Fig. [U The data for m ;d is well 
fitted by 



3.2 Sharpness of the transition 

For the term 'transition' to have any real meaning, it must 
be sharply defined in the N — ► oo limit, i.e. there exists 
a <j) c such that P = for all cj) < 4> c and P ^ for 
<p > <j) c . For finite N, however, this transition becomes 
'blurred' and the probability P[p^o] (0) that a system be- 
comes jammed (here identified by P ^ 0), varies smoothly 
from to 1 around </> m id ^ 4> c - Here we focus on the 
U co = 10~ 2 case, and as shown in Fig. [3] the expected 
trend of a sharpening transition as N increases is indeed 
observed. 

By fitting P[p^o] (0) t o the functional form 

P[p^o\{4>) = + tanh[(<£ - (j) m id)/o]} (2) 

it is possible to extract the width of the transition a and 
the point when p[p^o]((j)) = 4, m id, as a function of sys- 



0mid = 0oo + AN-" (3) 

with (jioo = C ([/ CO =1CT 2 ) = 0.2364 ±0.0017, A = 0.322 ± 
0.015 and v = 0.398 ± 0.018. This exponent v appears to 
be consistent with the corresponding repulsive-only value 
w 0.47 ±0.05 0. 

The width a is also a monotonic decreasing function 
of N, suggesting a sharp transition in the N — > oo limit. 
Although the statistics are poorer for these data and no 
convincing fit was possible, from the figure it is clearly 
power law decay or faster. We therefore feel justified in 
asserting that the transition is sharply defined in the infi- 
nite system limit. 

3.3 Coordination number and rigidity 

A cluster can be regarded as solid if its (free) energy in- 
creases whenever a non-trivial displacement field is ap- 
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plied to its particles, where 'non-trivial' means not rigid 
body translation or rotation of the cluster as a whole. If 
all contact forces are initially zero (as typically consid- 
ered in e.g. lattice models), the internal energy can only 
increase if the contact energies increase [3D]; for central 
force interactions as here, this corresponds to a change in 
the overlaps of contacting particles. The engineering term 
mechanism, meaning a displacement mode that does not 
alter the contact lengths, then coincides with the physics 
term floppy mode, meaning a zero-energy mode for which 
vanishing elastic moduli are expected. 

A simple counting argument, known as Maxwell count- 
ing, allows the number of mechanisms to be estimated 
from the network topology. As carefully detailed by Calla- 
dine [21], 

dN-N c + s = m+ ^d(d + 1) (4) 

with d the dimension, N the number of particles, A^ c the 
number of contacts, m > the number of mechanisms, 
and s > the number of linearly-independent stressed 
states. It is normal to express (j4]) in terms of the coor- 
dination number z = 2N C /N and to consider the limit 
N — > oo. Then, if we assume only one of m or s can be 
non-zero at the same time [see (iii) below], m > only 
when z < ^Maxwell = 2rf = 6 for d = 3, from which we infer 
the system is floppy. Conversely, s > when z > ^Maxwell 
(so there are 'redundant bonds') and we expect a finite 
restoring force for an arbitrary perturbation. 

Thus this counting scheme allows rigidity to be deter- 
mined from z alone. However, this analysis is only partial: 

(i) Transverse contact forces (i.e. friction), particle as- 
phericity, and non-pair (3 or more) particle potentials 
change the constraints per contact and degrees of freedom 
per particle. 

(ii) It refers only to the global-spanning rigid clus- 
ter; independent subsystems should not be including in 
the counting. For repulsive-only particles, such subsys- 
tems have been identified as individual rattlers 7,22]; for 
bond-diluted lattice models they are lattice animals that 
often occupy the bulk of the system [23 . 

(iii) Any form of linear dependency between the con- 
straint equations permits m > and s > simultane- 
ously, as demonstrated by Calladine [21]. Such stresses 
make some mechanisms non-floppy. 

(iv) The procedure tells us nothing about the stability 
of any solid state [TOlITT]. and thus may predict solidity 
for unstable clusters (which will spontaneously rearrange 
and hence are clearly not solid). 

(v) Ordered systems may respond anomalously to spe- 
cific strains. For instance, a rectangular lattice of unstressed 
Hookean springs is non-rigid to shears parallel to either 
of its lattice vectors, but rigid to all others [24] . 

(vi) Pairs of non-convex particles can make multiple 
contacts, modifying the counting procedure. 

(vii) There is a small error due to the rigid body trans- 
lation and rotation modes of the cluster as a whole, which 
must be subtracted from the degrees of freedom. For a free 
body this correction is \d(d +1). 



Plots of z a n = z versus 4> are given in Fig. [5] as well 
as the same quantity restricted to compressive and tensile 
bonds, ^compressive and ztensiie respectively, for various U co , 
including the repulsive-only case U co = 0. For U co = the 
Maxwell approximation is good, i. e. P becomes non-zero 
when z a n ss ^Maxwell- In fact equality can be reached once 
points (ii) and (vii) above have been corrected for [22] . 

The introduction of attractions does not appear to 
drastically alter the picture. For instance, the U co = 10 -2 
data, for which the transition occurs over a range 0.25 < 
<fi < 0.27 (see subsection 13.21 above), all three z's exhibit 
a knee around the transition range, which will presum- 
ably become a sharp discontinuity in the thermodynamic 
limit N — > oo. The knee corresponds to a sharp increase in 
•Ztensiic, confirming the transition to a tension-dominated 
state. 

Inspection of the raw data suggests z a u « 5.8 for jam- 
med configurations around <fi c with U co — 10 -2 , slightly 
less than the repulsive-only value w 5.86 [10 . Even taking 
into account the correction mentioned in (vii) (which is 
O(0. 01) here), this still leaves a significant deviation from 
^Maxwell = 6. The obvious explanation is that it is due to a 
small fraction of rattlers. However, for the attractive case 
here we must also consider the possibility of independent 
clusters of many particles. Indeed, the observation of z & \\ > 
for very low volume fractions with P = 0, which must be 
due to such clusters, suggests they can make a measurable 
contribution to z a n above <fi c as well. 

There is another possible contribution to the gap z < 
^Maxwell which cannot a priori be discounted. The exis- 
tence of self-stressed mechanisms mentioned in (iii) above 
permits rigid configurations for lower z than expected. For 
this to be significant, the particles must relax into a con- 
figuration in which some constraint equations are linearly 
dependent. This may appear to be a measure-zero event, 
but given the stabilising role of such states (see [21]), they 
may be attractors of the dynamics and hence arise with 
finite probabilility. For example, a string of attractive par- 
ticles pulled taught will naturally form colinear bonds, 
which are dependent constraints, and such a chain is rigid 
when under tension (but floppy when unstressed). 

4 Discussion 

The main finding of this paper has been the demonstra- 
tion of a well-defined transition to a low-density tensile 
solid in systems of athermal particles with attractive inter- 
actions. This extends the range of model 'jamming' tran- 
sitions that may shed light on the underlying nature of 
this process. Indeed, if the repulsive-only transition was 
to be crudely referred to as an athermal glass transition, 
the system here could equally crudely be referred to as 
an athermal gel transition. Just as in thermal gels (such 
as colloid-polymer mixtures, see e.g. |25p26j). the non- 
equilibrium solid forms at low densities, although it is not 
clear if there is any athermal equivalent of the glass-glass 
transition. 

Only one preparation procedure has been considered 
here, namely the simple quench from a random, T = oo 
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Fig. 5. The coordination number for all contacts 2 a u (upper 
panel, ZMaxwell = 6 also shown), and restricted to compressive 
and tensile contacts (lower two panels) for the given U co . For 
this N = 1000 system, the U co = 10 -2 transition occurs occurs 
over a range 0.25 < <fi < 0.27 (see Fig. [3}, around where the 
knee in the data occurs. For clarity, only the z&u = ^compressive 
data is given for the repulsive-only case U co = 0. 



state. Whereas repulsive-only particles produce a some- 
what homogeneous state, perhaps explaining their robust- 
ness to different preparation procedures and the near- 
ubiquity of 'random close packing,' attractive-particle ag- 
gregates are far more hetereogeneous, suggesting greater 
sensitivity to the means of preparation. Indeed, irreversible 
gelation of emulsions under shear has recently been ob- 
served [37], suggesting that jamming could be seen for 
even lower volume fractions than those observed here if 
the system is sheared during preparation. 

An original intention of this work was to determine 
the scaling behaviour of various properties, such as pres- 
sure or elastic moduli, near the transition point <p c , as it 
would be interesting to compare these to the dimension- 
independent rationals observed in the repulsive case [7j. 
In particular, to see if the visibly less compact nature 
of the network in Fig. [1] results in fractal spanning clus- 
ters at (j) c and irrational exponents. However, the use of 
standard minimisation algorithms (conjugate gradient and 
BFGS [17] ) failed due to the non-smoothness of the inter- 
action potential |T]), which has a discontinuous first deriva- 
tive at £ = (unlike the strictly repulsive case, which al- 
ways has a continuous first derivative, and for Hertzian 
contacts a continuous second derivative also). It proved 
impossible to make the two algorithms agree near to cf> c , 
and so the P data was discarded here (although P[p^o] {4>) 
seems to be robust). 

Future work could avoid problems with non-smooth 
potentials by employing a molecular dynamics algorithm, 
which would also allow the use of standard surface ten- 
sion interactions [16]. It would also be interesting to test 
the robustness with respect to particle polydispersity and 
system dimensionality, to see if the exponents remain d— 
independent. If the exponents did differ from the repulsive- 
only case, one possible interpretation is that the jammed 
configurations are indeed determined by elastic stability, 
since stability can be very different under tension than 



under compression(see e.g. [28]). In this context, it is par- 
ticularly interesting to note that rational, dimension-in- 
dependent exponents are completely typical of elastic in- 
stabilities. 
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